function [snr] = calcSNR(det, source, skyPos, psdType)

% Constants
MSOLAR_TIME = 4.92579497077314e-06;
PARSEC_SEC  = 1.0292712503e8;

Delta	= pi/2 - skyPos(:,1);
Alpha	= skyPos(:,2);
r	= source.r*PARSEC_SEC*1e6;

% Calculate the skyPos dependent factors
[at, bt] = CalculateAtBt(det, Alpha, Delta);

c    = cos(source.iota);
Phi0 = source.Phi0;
psi  = source.psi;

a1 = ( - (1+c^2)*cos(Phi0)*sin(2*psi)...
       - 2*c*cos(2*psi)*sin(Phi0) ) / r;

a2 = (   2*c*cos(Phi0)*sin(2*psi)...
       + (1+c^2)*cos(2*psi)*sin(Phi0) ) / r;
   
a3 = ( - (1+c^2)*cos(Phi0)*cos(2*psi)...
       + 2*c*sin(2*psi)*sin(Phi0) ) / r;

a4 = ( - 2*c*cos(Phi0)*cos(2*psi)...
       + (1+c^2)*sin(2*psi)*sin(Phi0) ) / r;


freq	= linspace(0.1, 3000, 10000);
df	= freq(2) - freq(1);
f0	= min(freq);
n0	= ceil(f0/df);
fLSO	= (1./sqrt(6))^3 / (pi*source.mTotal*MSOLAR_TIME);
nLSO	= ceil(fLSO/df);
Shf	= detnoisepsd(psdType, freq);
X	= dot( freq(n0:nLSO).^(-7/3), df./Shf(n0:nLSO) );

C1	= sqrt(5*source.eta/96) * (source.mTotal*MSOLAR_TIME)^(5/6) * pi^(-2/3);
C2	= sqrt((a1*bt + a3*at).^2 + (a2*at + a4*bt).^2 );
snr	= 2*C1*C2*sqrt(X);






